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Abstract 

In a series of recent papers (Adcock, Hansen and Poon, 2013, Appl. Comput. Harm. Anal. 
45(5);3132-3167), (Adcock, Gataric and Hansen, 2014, SIAM J. Imaging Sci. 7(3):1690-1723) and 
(Adcock, Hansen, Kutyniok and Ma, 2015, SIAM J. Math. Anal. 47(2):1196-1233), it was shown 
that one can optimally recover the wavelet coefficients of an unknown compactly supported func¬ 
tion from pointwise evaluations of its Fourier transform via the method of generalized sampling. 
While these papers focused on the optimality of generalized sampling in terms of its stability 
and error bounds, the current paper explains how this optimal method can be implemented to 
yield a computationally efficient algorithm. In particular, we show that generalized sampling has 
a computational complexity of O {M(N)\ogN) when recovering the first N boundary-corrected 
wavelet coefficients of an unknown compactly supported function from M{N) Fourier samples. 
Therefore, due to the linear correspondences between the number of samples M and number of 
coefficients N shown previously, generalized sampling offers a computationally optimal way of 
recovering wavelet coefficients from Fourier data. 


1 Introduction 

The motivation behind the papers 0 E E] is that, there are countless applications in which one 
is concerned with the recovery of a function from finitely many pointwise evaluations of its Fourier 
transform. To name a few examples, these applications range from medical imaging IMIMESl to 
electron microscopy [30] ? helium atom scattering |2E1 US] ? reflection seismology m and radar imaging 
m- More precisely, one is tasked with the recovery of an unknown function / G supported 

in a compact domain, from given samples of the form {f{io) : lj G H}, where Q is a countable set in 
and the Fourier transform is defined as 

ho = [ 

The set of sampling points Q either can be taken on an equidistant grid, which is known as uniform 
sampling, or can be an arbitrary countable set of certain density, which is known as nonuniform 
sampling. The listed papers described how, using generalized sampling, one can in principle optimally 
recover the (continuous) wavelet coefficients of a function, given its Fourier coefficients, and thereby 
directly exploit the superior approximation properties of wavelets. The present paper addresses the 
computational aspects of this generalized sampling method for both uniform and nonuniform sampling 
scenarios. 
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fCCA, Centre for Mathematical Sciences, University of Cambridge, UK (cmhsp2@cam.ac.uk) 
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1.1 The need for generalized sampling 

The classical approach to this recovery problem is to directly invert the Fourier sampling operator. In 
the case where consists of equispaced points, one can simply apply the Nyquist-Shannon sampling 
theorem to compute a discrete inverse Fourier transform on the given samples, while in the case where 
n consists of nonequispaced points, one can approximate / by inverting the frame operator associated 
with the Fourier frame : a; G [IHl HI]- Another approach to nonuniform sampling is so- 

called gridding m, a heuristic method which simply discretizes the inverse Fourier transform on a 
nonequidistant mesh by using certain density compensation factors. 

In practice, it is always the case that only a finite number of samples is available, i.e. card(r2) = 
M < oo, and hence, unless the original function is periodic and continuous, these direct Fourier 
approaches suffer from a number of drawbacks, including slow convergence rates and artifacts such as 
the Gibbs phenomenon. Therefore, these classical approaches are considered unsatisfactory (SUES]- 

In order to obtain a better reconstruction, one might seek an approximation represented in a 
different, more favorable system than the one used for sampling. In fact, a typical magnetic resonance 
(MR) image, for example, is known to be better represented by wavelets than by complex exponentials 
1371. The idea of basis change is not new in signal processing and it is typically referred to as 
generalized sampling. It can be tracked down to the work of Unser and Aldroubi [38L 125] . which 
was later extended by Eldar [MED]. Note that some of the early works of generalized sampling are 
sometimes referred to as consistent sampling. The idea of generalized sampling has also been applied 
in the work of Pruessmann et al. [21| and Sutton, Noll and Fessler [3S] for the recovery of voxel 
coefficients from Fourier samples and also in the work of Hrycak and Grochenig |26| for the recovery 
of polynomial coefficients from Fourier samples. Furthermore, the efficient recovery of trigonometric 
polynomial coefficients of bandlimited functions from its nonuniform samples, or equivalently, the 
recovery of Dirac coefficients of compactly supported functions from nonuniform Fourier samples, 
was considered in works by Feichtinger, Grochenig and Strohmer m- 

These approaches were formalized under the framework of Adcock and Hansen [5] in the setting 
of a Hilbert space "H, where it was shown that one can recover an unknown element / G "H in terms of 
any frame when given samples of the form ((/,'0j))jeN where is an arbitrary frame 

in li. More specifically, given the first M samples {(/, j = 1,..., A/}, one can approximate the 
first N reconstruction coefficients {(/, : j = 1 , ..., N} in a stable and quasi-optimal manner (see 

Section]^ for the definition) provided that M is sufficiently large with respect to N. 

The setting where forms a wavelet basis and forms a Fourier basis was considered 

in [H] and [5] in one and two dimensions respectively. It was shown that a linear scaling between the 
number of Fourier samples M and the number of wavelet coefficients N is sufficient for stable and 
quasi-optimal recovery. In other words, it suffices to take M = O (N) Fourier samples of a function 
in order to stably quasi-optimally recover its first N wavelet coefficients. A more general setting 
where forms a (weighted) Fourier frame was considered in [T], in one dimension, where it 

was shown that it is sufficient to take Fourier samples from the interval [—K,K] with K = O {N) in 
order to recover the first N wavelet coefficients. Therefore, as long as the number of samples M scales 
linearly with the sampling bandwidth K. the relationship M = O {N) is also sufficient for stable and 
quasi-optimal recovery in this more general setting. 

1.2 Contribution of this paper 

The method of generalized sampling can be recast as the linear system 
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which is to be solved for the least-squares solution {aj : j = 1,..., A^} that approximates the desired 
coefficients {(/, : j = 1,..., given the M samples {(/, '0^) : j = 1,. •., M]. In the general case, 

solving this system has a computational complexity of O {NM). In this paper we show that in the case 
of recovering wavelet coefficients from Fourier samples, the computational complexity of generalized 
sampling is only O (Mlog A^), since in this case the cost of applying the matrix from and its 
adjoint is only O(AflogA^). In fact, in uniform sampling as well as in nonuniform sampling where 
M = O (K), the computational complexity is simply O {NlogN), due to the linear correspondences 
derived in the aforementioned papers. Therefore, one can attain the superior reconstruction qualities 
via generalized sampling at a computational cost that is on the same order as the computational cost 
of the classical methods that are based on a simple discretization of the inverse Fourier transform. 

This paper will describe the computational issues relating to the recovery of boundary-corrected 
wavelet coefficients |15| when given either uniform or nonuniform Fourier samples. Although we shall 
only address the reconstruction of coefficients in one and two dimensions, the techniques described 
in this paper can readily be applied to higher dimensional cases. Generalized sampling may also 
be efficiently implemented with wavelets satisfying other boundary conditions (such as periodic or 
symmetric boundary conditions), however, here we consider the boundary-corrected wavelets of [15] 
since such wavelets preserve vanishing moments at the domain boundaries and form unconditional 
bases on function spaces of certain regularity on bounded domains. 

While in this paper we mainly focus on the linear recovery model (1.11, the computational aspects 
analyzed in this paper arise in various other nonlinear recovery schemes such as the ^i-minimization 
schemes introduced in [MlISlIl]. Namely, whenever one wants to recover wavelet coefficients from 
Fourier measurements, one needs fast computations involving the same matrix as the one appearing 
in (l.ll, as well as the fast computations involving its adjoint. Hence, the algorithms described in 
this paper can readily be applied to yield efficient implementations of these other nonlinear recovery 
schemes. 

We remark that although the efficient recovery of wavelet coefficients from Fourier samples was 
already considered in [24] . that work did not provide a comprehensive explanation of the imple¬ 
mentation using boundary-corrected wavelets and various types of Fourier sampling. Moreover, 
the purpose of this work is to provide a practical guide to the use of generalized sampling with 
wavelets and Fourier sampling. A MATLAB implementation of our algorithm is available at http: 
//www.damtp.cam.ac.uk/research/afha/code/. 


1.3 Outline 


Section describes the framework of generalized sampling. Section describes the special case of 
generalized sampling which produces wavelet reconstructions from Fourier samples. The wavelet 


reconstruction space and Fourier sampling space are defined in Section 3.1 and Section [3.2[ respectively, 
and additionally, in Section [3.3[ we recall the results providing the linear stable sampling rate for this 
particular pair of spaces. In Section]^ we deliver the detailed reconstruction algorithm for the one¬ 
dimensional case together with the analysis of its computational complexity. The two-dimensional case 
is presented in Section in a general setting, and the simplifications arising from uniform sampling 
are presented in Section]^ In Section we use the presented algorithm in a few numerical examples. 


1.4 Notation 

Let "H denote a Hilbert space with norm 11*11 and inner product (*, •). A set {0m}meN ^ ^ forms a 
frame for H if there exist constants A,B > 0 such that 

A\\ff<Y,\{f,^^)f<B\\ff, ( 1 . 2 ) 

mCN 

For any subspace C 'H, let Qj- denote the orthogonal projection onto T. Throughout the paper, 
we shall consider the Hilbert space of square-integrable C? functions which are compactly supported 
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on the cube [0, l]'^ C in the physical domain, for dimension d E 'N. For a point in the frequency 
domain cj G let 

e^{x) = ^ ^ 

Hence, for a function / G i3^([0,1]'^), we can write 

/(^) = 

for the Fourier transform of / at uj. In particular, if ( |1.2[ ) holds with fbat {eLij^}m,eN 

is a Fourier frame for i3^([0, For a matrix A G we use the notation A = {Am,n)n^m=i 

and denote its pseudoinverse by A^ and its adjoint by A*. 


2 Generalized sampling 

This section recaps the main features of the generalized sampling method described in laizi. Let 7Z,S 
be closed subspaces of H such that TZDS^ = {0} and 7l-\-S^ is a closed subspace of H. Let 
and {V'mImgN frames for TZ and S respectively. For A^, M G N, let 

IZn = span {ifri : n = 1,..., iV} , Sm = span : m = 1,..., AL} 


denote the finite dimensional reconstruction and the sampling space respectively. For an unknown 
function / G 'H, we seek an approximation f E IZn to f using only the finite set of samples 
{(/,'</’m) : rn= 

It was shown in [7] that for each iV G N, there exists mo E N such that for all M > mo, there 
exists a unique element / G 7Zn satisfying 

{SmI^^Pu) = {SMf,Pn), n = l,...,N. (2.1) 

where Sm • H Sm, 9 ^ The approximation / is referred to as the generalized 

sampling reconstruction. Furthermore, this reconstruction satisfies the sharp bounds 

||/-/||<C^V.M||Q7^./-/||, ll/l <C^.m||/||, (2.2) 


where 


uEI^N 

lkll=i 


u\\ > 0. 


Note that the approximation / achieves, up to a constant of Cn,m, the same error bound as Q^Zj^f, 
which is the orthogonal projection onto IZn und therefore the best possible approximation to / in 
TZn' Hence, / may be said to be quasi-optimal. 


Due to (2.11, by defining the generalized sampling operator 


N 


M 


GlN’M] : ^ a ^ an<Pn, i’, 


(2.3) 


n=l 


the solution / can be written as 


N 

1=1 


where a = G is the least squares solution to 


(«) = ((/, 


(2.4) 
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Note that this is the same as the system written in (1.1). Furthermore, the condition number of the 
operator which indicates reconstruction stability, can be shown to be 


Dn,m 


M 


inf 

U^IZ-N 
> m=l 


-1/2 

> 0 . 


The accuracy and stability of generalized sampling are therefore governed by Cn,m and Dn^m 
and it is imperative to determine the correct scaling between the number of reconstruction vectors 
N and the number of samples M to ensure boundedness of these two values. For a given number 
of reconstruction vectors N and some 0 G (0,cx)), the stable sampling rate is defined as a minimal 
number of sampling vectors M providing a stable and quasi-optimal reconstruction. Formally, the 
stable sampling rate is given as 


B{N,0) = min{Af G N : ma,x{CN,M, Dn,m} < 0} . 


(2.5) 


The stable sampling rate for particular choices of wavelet reconstruction spaces and Fourier sampling 
spaces is discussed further in Section [T3| 


2.1 The computational challenges 


There are two aspects to the implementation of generalized sampling in (2.1): 


1. We need to understand the stable sampling rate M = (3(N,0) required for quasi-optimal and 
stable reconstructions. 

2. Having d eterm ined the appropriate finite section matrix it remains to solve least-squares 

problem (2.4). Note that given any matrix G G there are many iterative schemes for 


computing the least-squares solution G^g. For example, one can apply the conjugate gradient 
method to the corresponding normal equations |36) . The efficiency of these iterative methods 
is always dependent on the computational complexity of applying G and G*, and in the worst 
case, this is O (NM). 

In the case of recovering wavelet coefficients from Fourier samples, the first challenge detailed above 
has been resolved in [H] where the sampling vectors form a basis, and in [J, where the sampling vectors 
form a (weighted) Fourier frame. The main results from these works will be summarised in Section 


3.3 


The purpose of this paper is to show that the computational complexity of applying 
its adjoint is 0{M\ogN). Thus, because of the linear correspondence between the 

number of samples and reconstruction coefficients revealed in [S] and [T], generalized sampling offers 
a computationally optimal way of recovering wavelet coefficients from Fourier data. 


Remark 2.1 (Generalized sampling and minimization) Based on ideas of compressed sens¬ 
ing introduced in [Mni], to recover j ^ C? from a reduced number of noisy Fourier measurements 
y = {f{ujj) '• j £ 0,} rf, where ||? 7||^2 < S and is some finite index set, it was proposed in 0 E] to 
solve 

X G arg inf ||-2:|L subject to \\Gqz — y|L < S, (2.6) 

where 

is a semi-infinite dimensional generalized sampling matrix. Solutions to this minimization problem 
were analysed in detail in [3133], and we refer to those works for theoretical error estimates. However, 
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we simply mention that if H = {j : j = —M,..., M} and (5 = 0, it was proved in [33] that if 


is a wavelet basis with sufficiently many vanishing moments, then any solution x to (2.6) will satisfy 


Ji jeN 


/-E 


Xjifj 


3 = 1 


= o 




L2 


.3=N' 


where N' = cN for some constant c G (0,1] which is independent of N. Thus, the reconstruction 
error is governed by the decay in the wavelet coefficients of /, rather than the decay in its Fourier 
coefficients. Computationally, one can show that any solution to 


inf llalLi subject to 




< S 

p 


(2.7) 


converges in an sense to solutions of (2.6) as ^ oo (see 0) and for the Fourier-wavelets case, 
it suffices to let K = O {M) [3]. There are again numerous algorithms such as EniEHj for solving 
problems of the form and, similarly as for ( |2.1[ ), the efficiency of these algorithms amounts to 
the efficiency of applying and its adjoint. 


3 The wavelet reconstruction from the Fourier samples 

As mentioned before, in this paper we consider the particular reconstruction space consisting of N 
boundary-corrected wavelets and the sampling space consisting of AI Fourier-type exponentials. In 
this section, we define the corresponding spaces and give the stable sampling rate for this pair. 

3.1 The reconstruction space 

We shall consider the Daubechies wavelets of a vanishing moments, on the interval [0,1], with the 
boundary wavelet construction introduced in m- Following convention, we assume that the scaling 
function (p and the wavelet ^ are supported on [—a -h l,a] and denote boundary scaling functions at 
the endpoints 0 and 1 by 

k = 0,...,a-l 

and boundary wavelets at the endpoints 0 and 1 by 

V’fc, V’L k = 0,...,a-l. 

We denote the boundary-corrected scaling functions on the interval [0,1] by 

( 2^/2^(2^a: - k) a < k < 2^ - a 
2-^/^(/)^(2-^x) 0 < A: < a 

[ 2^'/2^i,_^_^( 2^-^) 2-^'-a<A:<2E 

and similarly for the wavelet functions 

One dimension. Here we define the reconstruction space for dimension d = 1, consisting of 
boundary-corrected wavelets on the interval [0,1] of a vanishing moments. Let J be such that 
J > log2(2(^)- The set 
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forms an orthonormal basis for i3^([0,1]). Furthermore, by defining 


= span = 0,..., 2-^ - l| , = span | : A: = 0,..., 2-^ - 1 


we have that -C^([0,1]) = © ^©j<j We now order elements of (3.1 

of wavelet scale, and we denote by W^r the space spanned by the first N wave. 
R> J > log2(2a) and N = 2^ we have the reconstruction space 


Wn = V 


_ w[o,i] 


© W 


[ 0 , 1 ] 




N - yj ^ yy j 

R _T/[0,l] 


) in increasing order 
ets. Thus, for fixed 

(3.2) 


Note that, by construction, dim(>VAr) = N = 2^ and also Wn = Vr 
For a function / G Wat we can write 


R-l 2^-1 


and also 


/(®) = J2 (®) + Z! Z! 

fc=0 j=J k=0 


/(®) = Z CR,k^Rl'ix) 

fc =0 

for some scaling coefficients and some detail coefficients dj^k- We recall that, given the scaling 
coefficients {cR^k : A; = 0,..., 2^ — 1}, it is possible compute the scaling coefficients {cj^k • ^ = 
0,..., 2^ — 1} and detail coefficients {dj^k ■ k = 0,..., 2-^ — 1, j = J,..., — 1}, and vice versa. 
This can be done by the discrete boundary-corrected Forward Wavelet Transform (FWT), which 
we denote by W and by W in two dimensions. The reverse operation is performed by the discrete 
boundary-corrected Inverse Wavelet Transform (IWT), denoted by IF“^, in one, and by in two 
dimensions. 


Two dimensions. The boundary-corrected wavelet basis for £^([0,1]^) is defined as follows. We 
introduce the following two-dimensional functions: 


As before, take J > log2(2a), and let 

T/ = ■Q<ki,k2<2^ -1} 

and 


.[ 0 , 1 ] 


[ 0 , 1 ] 


,[ 0 , 1 ] 


'Y'l _ 


■0<ki,k2< V - l} , j e N, j > J, i = 1,2,3. 


The set 


TTjU I U = 1.2,3} 

3>J 


(3.3) 


forms an orthonormal basis for L^([0,1]^). We now order the basis elements of (3.3) in increasing 
order of wavelet scales so that we can write 


( TO Ti 


{V^ni,n2}ni,n2eN2 — 


^2 y3 


^2 

J+1 


n+i --A 
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Let denote the space spanned by the first N x N wavelets via this ordering, so that >VAr 2 = 
span{(^^^ ,12 : 0 < ni,n 2 < — 1}. For N = 2^ and R> J > log2(2a), we have 

Wats = span Tj 0 span T* j = span T% (3.4) 

which is the two-dimensional reconstruction space. 


3.2 The sampling space 

Let iiT > 0 be a sampling bandwidth and Z^: C be a closed, simply connected set in the frequency 
domain such that max^^^^ = K. For finite M = M(K), let 

Hm = {wm : m = 1,... ,M} C Zx, 


be a given sampling scheme, a set of distinct frequency points in Zk which are permitted to be 
completely nonuniform. The sampling space that we shall consider in this paper is given by 

8 m = span{y^e ■ ^rn e Om} , (3.5) 

with the weights defined as measures of the Voronoi regions associated to the sampling points, i.e, 

{x) dx, (3.6) 

Jzi^ 


where 


= {v e Zk - yn m,\v - < |^ - Wn|^i} • 

The use of Voronoi weights is common in both practice and nonuniform sampling theory [Mini- 
As shown in [2], it is crucial that the sampling scheme Qm satisfy an appropriate density condition 
in the frequency region Zk’, i-e. it is crucial that 


^Zk{^m) = sup inf \v -uj\^i < 1/2, 
v^Zk 


(3.7) 


but due to use of weights, Qm is allowed to cluster arbitrarily as K, M{K) oo. Namely, the density 
condition ensures that weighted exponentials in Sm give a weighted Fourier frame for £^([0,1]^) as 
K,AI(K) —^ 00 , i.e. (1.2) is satisfied with (j)m = y/firn^ujm • This, on the other hand, ensures stable 
sampling. 

The particular examples of nonuniform sampling schemes that we consider in the present paper 
are jittered and log sampling schemes in one dimension, and spiral and radial sampling schemes in 
two dimensions. The construction of these particular sampling schemes, such that they satisfy density 
condition (3.7), was described in [T] for the one-dimensional case, and in [2] for the two-dimensional 
case. 

A special case of sampling that we also consider in this paper is uniform sampling. Here, the 
sampling scheme is taken on an equidistant grid with a fixed spacing. For a given spacing e G (0,1], 
the uniform sampling space in one dimension is defined as 


8m = spanjecefc : k = - [M/2],..., [M/2] - 1} 


(3.8) 


and in two dimensions as 


£m^ = span{eie2e£ifcie£,fc, : - \M/2] < < \M/2] - 1}, (3.9) 

for some G (0,1], i = 1,2. It is known that exponentials in the uniform sampling space lead to a tight 
Fourier frame for i3^([0,1]^). As we shall see, for this special sampling scenario, the two-dimensional 
computations can be considerably simplified. 




Remark 3.1 Since the clustering of the sampling points is allowed, the stable sampling rate in the 
case of nonuniform sampling is defined by 


e{N,e) = min{_ft: > 0 : max {Cat,£’ jv,M(if)} < 0} , 


(3.10) 


for a given number of reconstruction vectors TV G N and some 0 G (0, cx)). Therefore, rather than a 
minimal number of sampling vectors M, as it is the case in (2.51, here we seek for a minimal sampling 
bandwidth K that provides a stable and accurate reconstruction. 


3.3 The stable sampling rate 

To emphasize the benefits of using generalized sampling to recover wavelet coefficients from Fourier 
samples, in this section, we recall the main results of mm and present a numerical example which 
demonstrates the superior convergence rates that generalized sampling offers. 

To begin, generalized sampling is particularly effective in the regime of recovering wavelet coeffi¬ 
cients from Fourier samples because a linear scaling between the number of Fourier samples M (or 
sampling bandwidth K) and the number of wavelets N is enough for stable and accurate recovery via 
generalized sampling. In one dimension we have the following. 


Theorem 3.1 ([SI H]). Let Wat be the reconstruction space consisting of N boundary-corrected 
Daubechies wavelets given by {3.2). 


(*) 


Let the sampling space £m be given by Then for any 7 G (0,1) there exists a constant 

Co = <^ 0 ( 7 ) such that 

min{M : inayi{CN,M^D^^M} < 7 } < cqN. 


(u) 


More generally, let LIm be contained in the interval [—K, K] and such that (5 = <5^ 
1/2, i.e. let the sampling space 8 m be given by (3.5). Then for any j ^ (0,1 — 5) 
constant cq = 09 ( 7 ) such that 


[-Ar,A'](^M) < 
there exists a 


min |i^ : max {Cat^m, ^ 


This result characterizes the Fourier samplings which permit stable and accurate recovery in spaces 
of boundary-corrected wavelets. Both parts of this theorem provide estimates on a stable sampling 
rate of the type 

e{N,c^) = o{N), 


where 0 is defined by (2.5) for uniform sampling, and by (3.10) for nonuniform sampling. Let us 


mention here that similar higher-dimensional results were obtained for the uniform Fourier samples 
in [S], while for nonuniform sampling, the higher-dimensional case remains an open problem. 

The result given in [32] states that if / G 1], s G (0, a), then ||/ — QvViv/ll — 0{N~^), where 
TL^ denotes the usual Sobolev space and Wa^ is the space of N boundary-corrected wavelets with a 
vanishing moments. This result in combination with Theorem |3.1| gives the following. 


Corollary 3.2 (O [I])- Let Wn be the reconstruction space consisting of N boundary-corrected 
Daubechies wavelets of a vanishing moments {3.2). Let f G "^^[0,1] with s G (0, a) be an arbitrary 
function. 


(i) Let the sampling space Sm be given by {3.8), then the generalized sampling solution f imple¬ 


mented with M = 0(iV, 7 ) samples satisfies 


\\f-f\\=0{M- 
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Let the sampling space Sm given by (3.5), then the generalized sampling solution f imple¬ 
mented with the sampling bandwidth K[M) = 0(A^, (1 + 5)/(l — 6 — 7 )) satisfies 


\\f-f\\=0{K-‘‘). 


Hence, up to constant factors, generalized sampling obtains optimal convergence rates in terms 
of the number of sampling points M (or the sampling bandwidth K), when reconstructing smooth 
functions with boundary-corrected Daubechies wavelets. 

These superior convergence rates given by Corollary |3.2| are depicted in Figure Namely, using 
an example of a continuous, nonperiodic function, we compare the convergence rates of generalized 
sampling with boundary-corrected Daubechies wavelets to the suboptimal convergence rates of the 
simple direct approaches based on the discretization of the Fourier integral. 


Equispaced sampling Jittered sampling Log sampling 



Figure 1: A nonperiodic continuous function f{x) = a:^cos(37ra:)x[o,i](a;) is reconstructed from pointwise 
samples of its Fourier transform taken on an equispaced grid with e = 1 (left), on a jittered scheme with jitter 
0.1 (middle) and on a log scheme (right), where the last two satisfy 5 < 0.97. Reconstruction is performed 
via GS using different types of boundary-corrected Daubechies wavelets: Haar, DB2 and DBS, and also 
via truncated Fourier series (TFS) in the uniform case, and gridding in the nonuniform cases, with density 
compensation factors defined as in (3.61. 


4 Computations in one dimension 


Let a be the number of vanishing moments of Daubechies wavelets, R > \o^ ^(2a ) the finest wavelet 
scale and N = 2^ the dimension of the reconstruction space Wn defined by {3.2|. Let M > N, and 
let Qm = {^m m = 1,..., M} be the given set of sampling points with associated Voronoi weights 
p^rn defined as in (3.6l, and let Sm be the sampling space defined by (3.51. For these choices of the 


reconstruction and the sampling spaces, generalized sampling reduces to the weighted least-squares 
system 

(4.1) 




M 


for a given vector of Fourier samples ((/, and solving for an unknown vector of boundary- 

corrected wavelet coefficients a = (an)^=i of a function /. 

Remark 4.1 As explained previously, efficient implementation of generalized sampling leans on the 
efficient implementation of the forward and adjoint operations, and (G^^’^^)*, which we 

describe in detail below. They can be implemented as a function handle in Matlab, for example, 
and then used in Matlab’s function Isqr for iterative solving of the least squares system (4.11. 
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From (2.31, it follows that, for spaces Wn and 6m^ given a G and ^ G C^, the forward 
operation can be written as 

V M 


N-1 


13 = (XI an‘Pn,y/Jh^e^J 


(4.2) 


and the adjoint operation can be written as 


n=0 


M 


m=l 


7 := (Gt^’""F(C) = (E 


m=l 


N-1 


n=0 


(4.3) 


We now describe how these operations can be computed efficiently by using the following operators: 

i) For the set of frequencies Qm and the corresponding set of Voronoi weights fim > 0, the diagonal 
weighting operator V = —>■ is given by 


It 


Vh) = (x4^7m)El . 7 e C". 

) For the set of freqnencies Qm, the operator F = is given by 

s M 


N-a-1 






k=< 


7 G C 


N 


(4.4) 


(4.5) 


m=l 


Hi) For the set of freqnencies Qm and the scaling function 0, the operator D = 
is given by 

/. \ \M 

(4.6) 


G C^. 


°«)=(K“^)c...):l. c 

For the weighting operator we have V* = V. The adjoint operator of F is F* : given by 

7 ^ J2m=i CmCo;^ k = a, — a—1 


(F*(C))fc = < ^ 

y 0 otherwise 

and the adjoint operator of D is D* : —)■ given by 

M 


C G 


- \ M 




(4.7) 


(4.8) 


Now we can analyse the operations (4.2) and (4.3). We first consider the forward operation. Given 
a G C"^, by the definition of discrete IWT, the equation (3 = a) is equivalent to 

iV-l 

Pm yJThn. ^ ^ dfc ipR,k ’ ^ 1? • • • i 

k=0 

where d = W~^{a) G C^. Since the Fourier transform of the internal scaling function can 

be written as 

(^) (-^) : k = a,...,N-a-l, 

by using definitions of operators F and F, we get 

Pm = -J=Pm = ^ E ( 9 ^) + {D (F(a)))^+ ^ E = M. 

VPm ViV ^^^k=N-a 
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Once /3 has been computed, it is left to apply the weighting operator and get (5 = V{P). 

For the adjoint operation, to compute 7 = {() given ( G we first apply the weighting 

operator and set = i^(C)- Then, similarly to the forward operation case, one can check that 
7 = VF “^7 and are related by the following equations. 


M 


(^)> k = 0,...,a-l, 7fc = ^ Yj U^N-k-i ('^)> k = N-a,...,N-l, 


m=l 


1 


M 


m=l 


and 


7fc 




k = a .— a — 1 . 


Vn 

m=l 

Note that, by using adjoint operators D* and F*, this last part can be written as 


7. = (f(D-(0))^. k 


= a,..., iV — a — 1. 


These computational steps, which we summarize below, lead to the efficient algorithm for forward 
and adjoint operations, and therefore to to the efficient algorithm for solving the weighted least-squares 
system (4.1 1 . 


The one-dimensional algorithm 


Precompute the Voronoi weights (fiTri)m=i and pointwise measurements of the Fourier transforms of 
the three scaling functions: 


HWL- mrL- mYL- . 

The forward operation 

GivenaeC^,= a) can be obtained by applying the following steps. 

1. Compute the scaling coefficients d = W~^{a), where W~^ is the one-dimensional discrete 
boundary-corrected IWT. 

2. Compute contributions from the boundary scaling functions: 

M 


M 


M 




1 

Vn 


k=0 


/3' = 


m=l 


1 

yiv 


E a.*-.-. Cw) 


k=N- 


M 


m=l 


3. Compute contribution from the internal scaling functions: 


3.1. Apply F to d to get a = F (d), where F is defined by (4.51. 

3.2. Apply F to d to get = F(d), where D is defined by (4.61. 

4. Compute 


5. Apply V to compute P = V{P)^ where V is defined by (4.4l. 
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The adjoint operation 

Given ( G C^, 7 = (C) can be computed as follows. 


1. Apply the weighting operator V and set = V{Q. 

2. Compute the coefficients of the boundary scaling functions: 


7fc = 


Vn 


M _ M _ 

^ ) , jN-k-1 = ^ ^ Cm0fc ) 

m=l ^ m=l 


for k = 0,... ,a — 1. 

3. Compute the coefficients of the internal scaling functions: 


3.1. Compute = D*{Q, where D* is defined by (4.8 1 . 


3.2. Compute jk = [F* 


^ k = a — — a — 1, where F* is defined by (4.7). 


4. Compute 7 = 1 ^( 7 ), where VF is a discrete one-dimensional boundary-corrected FWT. 


Remark 4.2 Note that when the sampling points are uniform with spacing e, the weights are simply 
fijn = € foi' all ^5 and the precomputation of Voronoi weights is not required. 

Remark 4.3 The above algorithm requires the precomputation of pointwise evaluations of the Fourier 
transform of the internal and boundary scaling functions. Note that for the internal scaling function 
(/), we may use the approximation 


N 

JJ mo(2“-^0 ^ N ^ 00 

j=i 

where uiq is a trigonometric polynomial cni. A similar approximation may be used in the case of the 
boundary scaling functions. This is outlined in detail in the appendix. 

Remark 4.4 On evaluating the reconstructed signal: Recall that in solving 

for an appropriate M = O (N), we obtain a which is an approximation of the first N wavelet coeffi¬ 
cients of / and the reconstructed sigual is f = 'l2f=i To evaluate / on the grid points {k2~'^)‘l^-^ 

for J G N, it suffices to evaluate each (j) on these grid points and we may do so by either implementing 
the cascade algorithm m or the dyadic dilation algorithm m- 

Computational cost of the one-dimensional algorithm. Let us analyse the computational 
cost of the forward operation. The adjoint operation can be analysed similarly leading to the same 
computational cost. The computational cost of applying the discrete boundary-corrected IWT is 
0{N). The cost of step 2, involving boundary scaling functions, is O(aM). For step 3.1, the key 
point is to observe that F is simply a restricted and shifted version of the discrete (nonuniform) 
Fourier transform, and thus its fast implementation (FFT/NUFFT) can be used when computing 
F(d). Hence, the the cost of step 3.1 is O {Mlog{N/e)) in the uniform case and O (Flog(iV) -h JM) 
in the nonuniform case, where L is the length of underlying interpolating FFT for NUFFT, and J 
is the number of interpolating coefficients (typically J = 7) [22]. Finally, the cost of the diagonal 
operations in both steps 3.2 and 5 is G (M). Therefore, given that e ~ 1, J ~ a and L ~ M, the total 
cost is essentially O (aM -j- Mlog(A^)). 
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5 Computations in two dimensions 

Again, let a be the number of vanishing moments, R > log2(2a) the finest wavelet scale and N = 2^, 
Let WAr 2 be the reconstruction space defined by (3.41. Let M > N, and let 


Om = { 






be the set of sampling points in R^, which we write as Qm = ^m) correspondingly. Let Sm be 

the associated sampling space defined by (3.51, with the Voronoi weights defined as in (3.6 1 . In 


this case, the least-squares system (2.4) becomes 


N 


M 


( /^m( 'y ^ 0^ni,n2^ni,n2 5 j ~ i ^ujr^})jn=l ’ 

V ni,n2 = l /m=l 

If we apply the two-dimensional boundary-corrected IWT to the matrix of wavelet coefficients a G 
so that d = G we get 


N 


M 


)^2 = 1 


m=l 


Since we can write the following algorithm. 


.[ 0 , 1 ] 


The two-dimensional algorithm 

Precompute the vectors and 



The forward operation 

Given a G I3 = G can be obtained by applying the following steps. 

1. Compute the scaling coefficients d = W^^(a), where is the discrete two-dimensional 

boundary corrected IWT. 

2. Compute contributions from the corners (the boundary scaling functions in the both axes): 


^00 

^01 

/310 


a —1 a—1 




fcl=0 fc2 = 0 
a-l N-1 


N 


M 


m=l 


E E 

fci=0 k2=N—a 
N-1 a-l 

E E <afcl,/c20}v-fci-l 

ki=N—a ^ 2=0 


m \ 11 

^iV-fc2-l 




N 

/ ,2 

m \ 20 I 


N 


4>l 


N 


M 


m=l 

M 


771 = 1 


A ^-1 A ^-1 / 1 \ / 2 \ \ 

ki=N-ak 2 ^N-a ^ ^ ^ / m=l 
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3. Compute contributions from the edges (the boundary scaling functions in only one of the axes): 


a—1 


v'’—a 

- a— 1 


fci=0 

1 

= ^ y .1 

fci=iV-a 


fc 2=0 

a— 1 


^intl _ — \ Dq^i (hD^'2 fki Fq\ (d. ^ 2 ) 




where F and D are defined by (4.51 and (4.6), respectively. 

4. Compute the contribution from the internal scaling functions: 

4.1. d = Fom (a), where Fn^ ■ ^ is such that for each 7 G 

M 


/ ^ N-a-1 . 

Fn„(7) = I F ( ' 

\ ki,k2 = a ^ 


{ki,k2) 

N 


4 . 2 . = D^.,Dn^Ja). 


5. Compute 13 = 13°° + 13°^ + I3^° + + /30“* + /3i“* + /3““ + /3“i + /3 


m=l 


jintint 


6 . Apply V to get /3 = V(I3), where V is defined by (4.4). 

The adjoint operation 

Given C G 7 = (0 G C^’'^ can be computed as follows. 

1. Apply the weighting operator V and set ^ = V{Q. 

2. Compute the scaling coefficients at the corners 


M 


IkiM — 'y^ Cini'k-i 


m=l 

M 


N 


m 1^0 

kr 


N 


lki,N-a-\-k2 — 


m=l 

M 


UJ: 


N 




a—k2 — 1 


N 


lN-a-\-ki,k2 — 'y Cm$l_ki-1 TV 


1 


M 


7N-a-\-ki,N-a-\-k2 


TV y, ^^^a-ki-1 jY 
m=l 




'a — k^ — i 


N 


for /ci, /c 2 = 0 ,..., a — 1 . 
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3. Compute the scaling coefficients at the edges 


7^1 ,^2 

1 

Vn 

((f„ 

■J 



J 


^ki,k2 

1 

\/]V 

((f„ 




-fc2- 

>'Gr 

* 

— a — 1, 

k2 = ' 

0, . . . : 

,a- 

1 and 




7fcl,fc2 ~ 

1 

((f„ 


(Dn^ 

\ “m 




7fcl,fc2 ~ 

1 

((f„ 

■j 

\ “M 

<Pa-ki-l J \ 


- 1, h = 

--a,.. 

.,N- 

- a - 

- 1, where and D* 

are defined 


respectively. 

4. Compute the scaling coefficients of the internal wavelets 

4.1. (0- 

4-2. 7fci,fc2 = (^F* ( 4 , 0 ))^ ^ 


= a — 1. N — a — 1. 


5. Compute 7 = W( 7 ), where W is the discrete two-dimensional boundary corrected FWT. 


Computational cost of the two-dimensional algorithm. For the forward operation, the com¬ 
putational cost of step 1 is O {N^) and that of step 2 is (9 Step 3 is (9 {a{M -h Mlog(iV/e))) 

in the uniform case and O (a(M + Llog{N) -h JM)) computations in the nonuniform case. The cost 
of step 4.1 is basically the cost of the two-dimensional NUFFT or FFT, i.e, O log(iV^) -h J‘^M) or 
O (Af log(A^^/e^)). The cost of step 4.2 as well as step 6 is (9 (M). Hence, if we assume 1, J Qj 
and Af, the total cost is O (a^Af -h A/logA^^). The same cost holds for the adjoint operation. 


6 Simplifications in the uniform two-dimensional case 

Let us assume the setting from the previous section, but now the sampling points are given uniformly 
i.e, the sampling space Sm"^ is given by (3.91 where, for convenience and without loss of generality, we 
assume ei = €2 = 1. In this case, we have 


q[N'^,M^] .^NxN ^ I^MxM ^ aH.f( an:,,n^(Pni,n^,emiem^)\ 

\ ni,712 = 1 / 


{M/2']<mi,m2<\M/2'\-l 


Let us analyze the forward operation, i.e. the computation of ^ given a. As 

usual, we can apply the discrete wavelet transform to get 


A^-l 

/^7rii,T7i2 ( 'y ^ R,{ki,k2)i ^ 1 ^ 1 ^ 1712 ) ■) [-44/2] ^ Im-i | , |77l21 ^ A/ 

fcl,fc2=0 

where a = W^^a. Now denote steps 2-3 of the forward operation in the one-dimensional algorithm 
applied to the scaling coefficients by . Given the special structure of both the 
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reconstruction and the sampling vectors, i.e, by using the definition of ^R^(ki,k 2 ) we 

obtain 


N-l N-l 


N-l 


^mi,m2 ( E(E (E 1 ) — Vmi, 

fci=0 A:2=0 fci=0 


m 2 


where for each ki = 0,..., iV — 1 


_/v—1 \ r-^/21—1 

rM/ 21-1 _ / / . .[0,1] 




k2 =0 




m2 = — [1^4/21 


and for each m 2 = — \M/2],... ^\M/2] — 1 


/ ^{M/2'\-i _ 

\Vmi,m2)mi=-[M/2] ~ [[Iki ,m2 ) ki=0 J ' 


Hence we can write the following algorithm. 


The two-dimensional algorithm in the uniform case 

Precompute the vectors 


/ " /m\\r^/ 2 l-i f 70 f 

y\N))m=-[M/2-\ ’ \N))m=-\M/2^ 



\M/2'\-l 
m=— fM/2] 


/c = 1 ,..., a. 


The forward operation 

Given a G /3 = can be obtained as follows. 

1. Compute the scaling coefficients d = 

2. For each fci = 0,..., — 1, apply G^p to the row of d to obtain 


.N-l 


Let 


Pki — G(j) (^(d/i;^^/i;2)fc2=0 

/ /3J \ 


-M 


/3 = 


>NxM 


WnJ 


3. For each m 2 = — [M/2],..., [M/2] — 1, apply Gp to the m^^ column of 13 to obtain 

N-l 


0m2 — Gr 


Let 


V /fci=0 
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The adjoint operation 

Given C e 7 = 




(Q can be computed as follows. 


1. For each k 2 = — [M/2],, [M/2] — 1, let the column of C, indexed by k 2 be denoted by so 
C'"" = Apply lo every column of ( to obtain 


7= ((G0)*(rl"/"1) ••• (G^)*(Cl"/"l-^)) ec 


NxM 


2. For ni = 0,..., iV — 1, let the row of 7 indexed by ni be denoted by . Apply to each 

to obtain 


1 = 


gC 


NxN 




3. 7= w [^). 


Computational cost of the simplified two-dimensional algorithm. Again, for the first step 
of the forward operation we have the computational cost O (A"^) • Since the cost of each application of 
operation Gfj) is O {aM + Mlog(A/e)), and the cost of each transposing of an M-dimensional vector 
is O (M), we have that the cost of step 2 is O (A(aM + Mlog(A/e) + M)) and the cost of step 3 is 
O (M(aM + Mlog(A/e) + M)). Therefore, the total cost amounts to O (aM^ + M^ log(A/e)). Note 
that the cost of the non-simplified two-dimensional algorithm in this case is O -h log(A/e)^). 


7 Numerical examples 

Finally, we present a few numerical examples using the algorithms from Sections [4][^ MATLAB 
scripts for reproducing all the examples in this section are available at http://www.daintp.cam.ac. 
uk/research/afha/code/ 

Example 1. In Figure we reconstruct a one-dimensional continuous, but nonperiodic, function 

f{x) = — exp(a:cos(47ra:))cos(77rx)x[o,i](a:) -h sin(37ra:)x[o,i](x) (7.1) 

from its Fourier samples taken on three different sampling schemes: equispaced, jittered and log. We 
use GS with boundary-corrected Daubechies wavelets of order 4 (DB4), and compare its performance 
to the direct approaches based on a discretization of the Fourier integral (truncated Fourier series and 
gridding). While the direct approaches are apparently plagued by Gibbs artifacts and the reconstruc¬ 
tion quality evidently depends on the underlying sampling scheme, GS performs equally well for all 
three sampling schemes producing a numerical error of order by using only 64 reconstruction 
functions. 
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Equispaced sampling Jittered sampling 


Log sampling 





bO 

.B 

bO 

CO 





Figure 2: Function is reconstructed from its Fourier samples taken from the interval [—64, 64] according 
to three different sampling schemes: equispaced (with spacing e = 1, in total 128 points), jittered (with spacing 
e = 0.77 and jitter 0.1, in total 168 points) and log (with density S = 0.97, in total 653 points). In the upper 
panels, reconstruction is done via GS with 64 DB4 wavelets, while in the lower panels via truncated Fourier 
series (TFS) in the uniform case, and in the nonuniform cases via gridding with density compensation factors 
l |3.6[ l. The £2 error is written below each reconstruction. 


Example 2. Next, we reconstruct a two-dimensional function shown in Figure which is again 
continuous but nonperiodic. We use an equispaced sampling scheme in Figurej^ and a radial sampling 
scheme that satisfies density condition (3.7) in Figure]^ By this example, we demonstrate robustness 
of GS when white Gaussian noise is added to the Fourier samples. We also demonstrate how one 
obtains improved reconstructions with increasing wavelet order. 


Example 3. We report numerical error of reconstructions from Example 4 in Table We see how 
error improves as the wavelet order increases. As evident, the main advantage of the GS approach 
is the possibility of changing the reconstruction space and using higher order wavelets for better 
performance. 


Example 4. In Figure we demonstrate reconstruction of a discontinuous function via GS with 
DB4, and compare it to the direct approach, when using uniform samples. 


Example 5. Figure illustrates the advantage of combining generalized sampling with mini¬ 
mization from (2.71, as opposed to simply employing standard finite dimensional compressed sensing 
techniques where one assumes compressibility with respect to some wavelet bases and one solves 


X G arg min llzILi subject to 
zec^ 




(7.2) 
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Figure 3: f{x,y) = sin(57ra;) cos(37ry)x[o,i]2 y) 


TFS GS with Haar GS with DB2 GS with DB3 



Figure 4: Top-left corner close-ups of the reconstructed function of / in Figure [fusing uniform samples. 


11/-/II 

SNR 

TFS/Gr 

GS with Haar 

GS with DB2 

GS with DB3 

Uniform 

0 

7.71 X 10-^ 

4.13 X 10-^ 

3.71 X 10-^ 

8.11 X 10-^ 

30 

7.88 X 10-^ 

4.23 X 10^'' 

9.14 X 10^^ 

8.74 X 10^^ 

Radial 

0 

1.98 X 10-^ 

4.13 X 10-^ 

3.74 X 10-^ 

7.95 X 10-^ 

30 

2.93 X 10^^ 

4.28 X 10^'' 

1.07 X lO^"' 

1.08 X 10^^ 


Table 1: The C 2 error for the reconstruction of function /(a:, y) = sin(57rx) cos(37ry)x[o,i]2 y) via truncated 
Fourier series (TFS) or gridding (Gr), and via GS with 64 x 64 Haar, DB2 or DB3 boundary-corrected wavelets 
from noiseless (SNR=0) and noisy (SNR=30) pointwise Fou rier measurements taken on a uniform or radial 
sampling scheme from the region [—64, 64]^. See Figures 3|4 and[^ 
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Figure 5: 


Top-left corner close-ups of the reconstructed function of / in Figure using radial samples. 



Figure 6: A discontinuous function is reconstructed by truncated Fourier series and by GS using 256 x 256 
DB4 boundary-corrected wavelets from 512 x 512 pointwise Fourier measurements taken from on an equispaced 
grid. 
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for given (noisy) Fourier samples y. The reconstruction is then taken to be U'^^x. Here, Udf is the 
discrete Fourier transform, Udw is the discrete wavelet transform and Pq is a projection matrix setting 
entries which are not indexed by Q to zero. However, as explained in [5], any solution to this will be 
limited by the error of the truncated Fourier representation of /, whereas the error of the solution to 


(2.71 can be controlled by the decay of the underlying wavelet coefficients. 


In this example, we consider the recovery of the smooth function 

f(x) = cos(3a:) sin(5y) exp(-x - y)A'[o,i ]2 


given access to samples of the form P = /C 2 ) : ki, k 2 = —512,..., 512}. Suppose that we observe 

only 4.25% the samples in namely, those restricted to the star-shaped domain H shown on the top 
right of Figure Note that the reconstructions obtained by solving ( |2.7| }, with K = 256 and ( |7.2| ) 
are both given the same input samples, but there is a substantial difference in reconstruction quality 


due to the samples mismatch introduced by the use of the finite dimensional matrices in (7.2). Note 


also that, when subsampling from the first M samples of the lowest frequencies, the computational 
complexity of computing (7.2) is O (Mlog(M)) since the computational complexity of applying Udf is 
O (Mlog(M)), and the computational complexity of applying Udw is O (M). Thus, the computational 


complexity of computing (2.7) is no worse than the finite dimensional approach of (7.2). 


A Computation of the Fourier transform of boundary scaling 
functions 

For the precomputational step in our algorithms, we need to compute the Fourier transform of bound¬ 
ary scaling functions. It is well known that the Fourier transform of scaling functions can be approx¬ 
imated using the low pass filter of the wavelet system. The key ideas on how to do this are recalled 
in this section. First observe that the internal scaling functions satisfy the following equation: 

(j){x) = y/ 2 'y^^ hn4>{2x — n) 


which gives rise to the equation 


^{ 0 = 

i=i 


where mo(0 = This provides a recursive way of approximating cf)^ since (/!>(^) ~ 

n^i large N. The same principle can be applied to compute the Fourier transform 

of the boundary scaling functions, for completeness, we will write out the computations in the case of 
the left boundary functions. 

Recall from ca that the left boundary scaling functions satisfy the following equations for each 
A: = 0,..., a — 1 

( a —1 a+2fc \ 

^=0 m=a / 

where|HQ ,| andj/ig ^} are the filter coefficients associated with the left boundary scaling functions 
(these filters are available from http: //www.pacm.princeton. edu/~ingrid/publications/54. txt). 
So, 




-27rim^/2 
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Original 


Sampling mask 



Figure 7: The reconstruction of a smooth function (top left) from 4.25% of its first 512 x 512 Fourier coeffi¬ 
cients using a standard compressed sensing (CS) approach with DFT and DWT, and using the generalized 
sampling (GS) approach. The samples taken for both reconstructions are those supported on the start shape 
sampling mask (top right). The bottom row shows a zoom-in of the top left reconstruction. The error from 
the standard CS approach is 1.6 x 10~^, while the error from the GS approach is 4.7 x 10~^. 
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Let 


u= 

( Ko 

■ 

Hli ■ 

■ Ka-l \ 

, 

f 

hi,a 

0 0 ... . 

^l,a+l ^l,a+2 0 

• • O O 


K^a-1,0 





'^a-l,a + l . 

• • ^a-l,3a-2/ 


and 


viiO = 

(\ 

II 

(M 

/ g-27ra^ \ 

g-27r{a+l)^ 




^g27rz{3a-2)^^ 


It follows that for each j G N 


1^1 (0 






We remark also that since 0(0) = 1, f2(0) is the vector whose entries are all ones and the values vi(0) 
can be obtained by solving 

i;i(0) = [/(t;i(0)) + l/(z;2(0)). 


Note that (jfi- is smooth and so, for each ^ G M, ^ r;i(0) as j —>• oo. So, 


^^ 1(0 - (^1 m+J2^‘v (^)) 


as j —>• oo and for j sufficiently large, we may approximate vi{^) by 

j-i 


viiO^U^v, {0)) + Y^U^vU (^ 
/=0 ^ 


(A.l) 


(A.2) 


Thus, for a given G M, to approximate we first solve (A.l I to obtain r'i(O), and then we 

compute (A.21 for j sufficiently large. 

A similar approximation can be derived for the Fourier transforms of the scaling functions on the 
right boundary. 
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